Computing the Frequency Response of Biochemical Networks: A Python module

In this paper, a set of Python methods is described that can be used to compute the frequency response of an arbitrary biochemical network given any input and output. Models can be provided in standard SBML or Antimony format. The code takes into account any conserved moieties so that this software can be used to also study signaling networks where moiety cycles are common. A utility method is also provided to make it easy to plot standard Bode plots from the generated results. The code also takes into account the possibility that the phase shift could exceed 180 degrees which can result in ugly discontinuities in the Bode plot. In the paper, some of the theory behind the method is provided as well as some commentary on the code and several illustrative examples to show the code in operation. Illustrative examples include linear reaction chains of varying lengths and the effect of negative feedback on the frequency response. Software License: MIT Open Source Availability: The code is available from https://github.com/sys-bio/frequencyResponse.


Introduction
Systems biology is a discipline that is concerned with understanding the behavior of biochemical systems in terms of their parts.This includes how perturbations in external factors or how variation in protein expression levels influence system behavior.More broadly it is also concerned with understanding the operational principles and logic by which cells manage resource allocation and growth.In many respects, systems biology is similar to physiology but whereas physiology has both quantitative and qualitative aspects, systems biology is firmly quantitative in its approach.
The easiest question to answer in systems biology is how components contribute to the overall behavior of a system.For example, consider the simple biochemical pathway: which has three reaction steps catalyzed by three enzymes e 1 , e 2 , and e 3 , a boundary species X o which we assume to be fixed so that the system can reach a steady-state (also called the fixed point or equilibrium point in other disciplines), and two species, x 1 and x 2 which can evolve in time.If the rates of the three reactions are given by v 1 , v 2 , and v 3 respectively, using mass-conservation we derive the differential equations for this system as: (2) arXiv:2406.11140v2[q-bio.MN] 27 Jul 2024 Note there is no differential equation for X o because X o is a boundary species that is considered an input to the system.
As an input, we can in principle apply any kind of input signal to X o .The most common signal, however, is a step input, which essentially makes the input a constant.X o isn't the only possible input, the enzymes along the pathway can also be designated inputs since we can imagine up-regulating their expression.Theoretically, any parameter of the system can be designated an input, though in practice not all parameters can be manipulated experimentally.
In this article, the words parameter and input will be considered synonymous; sometimes I will use the word parameter instead of input because it will be more appropriate in those instances.
There are some obvious questions we can ask, for example, how does enzyme e 2 influence the steady-state concentrations of species x 1 and x 2 ?If we were to increase e 2 does that increase or decrease the steady-state level of x 1 ?In addition, what determines the degree of influence e 2 might have of x 1 and x 2 ?
One way to answer such questions is to build quantitative models using differential equations from which we can run computer simulations to obtain time courses.However, trying to understand a system this way is not easy unless the system is very simple, and even then subtleties can be missed.For more complex systems, understanding the system's behavior by relying only on computer simulation is surprisingly difficult.Instead, various theoretical and numerical tools have been developed to assist a modeler in unraveling why a system operates the way it does.
One of the most common methods is to focus on the steady-state and study perturbations around the steady-state.
There is however one issue that limits what we can do.Most models found in systems biology, other than the simplest, are nonlinear.For example, reaction rates might be governed by nonlinear equations such as Hill equations or other enzymatic rate laws such as Michaelis-Menten rate laws.In addition, many reactions, particularly in signaling networks, are bimolecular which can also lead to nonlinearities.Nonlinearity makes such systems mathematically intractable.
As a result, it is common practice to linearize1 such systems around the steady-state.This makes the mathematical analysis tractable however it does limit the analysis to regions close to the steady-state.Even with these restrictions, a wealth of knowledge can still be obtained from linearized models.This is the basis of theoretical approaches such as metabolic control analysis and the closely related biochemical systems theory (Savageau, 1972;Kacser, 1973;Heinrich and Rapoport, 1974;Sauro, 2018) To study nonlinear behavior other techniques are available such as the use of phase portraits, bifurcation analysis, and a host of methods if the system exhibits chaotic behavior.In general, however, understanding nonlinear systems is difficult.
In this article, the focus will be on one specific method, the frequency response, that can be employed when studying linearized models.
I should say a few words on some specific terminology, X o , from a control theory perspective, is called an input.In systems biology, such inputs are also called boundary species (Hucka et al., 2003).Species that are not classed as inputs, such as x in scheme (1), are called state variables in control theory or floating species in systems biology.Floating species will evolve in time as a result of the operation of the model.In a later section, I will use the terminology 'floating species' rather than 'state variable' and one should be aware of their meaning.

Frequency Response
Because the frequency response of a system is not well known by the systems biology community and surprisingly, by the dynamical systems community, it is worth explaining this technique in a little detail.
A common technique that is used in designing and understanding man-made devices is the frequency response.This technique has not found much use in systems biology (Mettetal et al., 2008) but is widely used in engineering fields.
Operationally obtaining the frequency response involves applying a sine wave signal at a chosen input and examining how that time-varying signal is manifested at the system outputs.For example, in the pathway shown in (1), we might apply a sine wave to the boundary species X o and observe how this sine wave propagates to the internal species x 1 and x 2 .Intuition suggests that both x 1 and x 2 will also begin to vary in response to the input and this is indeed the case.However, the components in the system will tend to distort the sine wave and the oscillations we see in x 1 and x 2 will not exactly match the input sine wave applied to X o .Of particular interest is that a linear model (or a nonlinear model that has been linearized), behaves in a special way when a sinusoidal signal is input to the system.In this situation, the system will distort the input sine wave by only changing the amplitude and phase of the sine wave.The amount that the system changes the amplitude and phase is dependent on the frequency of the input wave.Such changes can give useful information on the properties of the system.It is important to note that the frequency component of a sinusoidal signal is unchanged in a linear system.Once it has we can see that its amplitude is reduced and its phase is shifted.That is the peaks don't match the peaks in the blue input sine wave.The observed phase shift is due to the time it takes to fill and empty the intermediate pool x.
Computing the frequency response involves computing the steady-state of the system and applying a sine wave at a given frequency and amplitude at one of the inputs.Since the input is now varying, the system responds by following the variation, however, due to delays in the system and the potential for attenuation or amplification of the sine wave, there will be changes in the amplitude and phase compared to the sinusoidal at the input.By varying the frequency of the input sine wave we can observe how the system responds across a range of frequencies.
In practice, we don't physically inject a sine wave into a computational model but can instead be done mathematically.For a physical system such as an electronic circuit, the frequency response is often obtained by actually injecting a sine wave into the circuit.
The mathematical techniques for obtaining the frequency response are well-known in engineering fields such as electrical and mechanical engineering and it is possible to those these techniques directly.However, there is an advantage to modifying these slightly by exploiting specific biological information such as the stoichiometry matrix and the reaction step elasticities (Kacser, 1973;Sauro, 2018).Ingalls (Ingalls, 2004) was the first to examine the frequency response from this perspective.One advantage of this approach is that it is possible to implement algorithms for computing the frequency response for any kind of biochemical network.
In this article, I will describe a Python package that allows a user to investigate a model, either derived from an SBML (Hucka et al., 2003) file or an antimony description (Smith et al., 2009) and compute the frequency response.Moreover, users can select any input and observe the frequency response on any output.This allows the user to obtain all combinations of inputs and outputs should they need to.The package uses libroadrunner (Welsh et al., 2023) as the simulation engine.In principle, however, any equivalent simulation engine will work so long as it can compute the steady-state, the elasticities, the reduced stoichiometry matrix, and the link matrix (Reder, 1988).

Theory
The frequency response of a system is based on transforming the system into the frequency domain.This assumes that the model is cast as a set of ordinary differential equations.The following description is based on the work by Ingalls (Ingalls, 2004).An additional publication was also published shortly afterward by Rao et al (Rao, Sauro and Arkin, 2004) who described a similar approach.The tutorial by Schulthess et al Schulthess et al. (2018), has a useful flowchart in their Figure 1 that describes some elements of the workflow.
For convenience, the symbols and dimensions used in the discussion are shown in Table 1.The starting point for the analysis is the system equation (Heinrich and Schuster, 2012;Hofmeyr, 2001) In this equation, N is the n × r stoichiometry matrix containing the stoichiometric coefficients of the floating species, with n being the number of floating species, and r the number of reactions.v is the r dimensional vector of reaction rates; x is the n dimensional vector of floating species and p the p dimensional vector of inputs.The p vector will depend on the application.All these are described in Table 1.
In metabolic control analysis, the inputs are often the vector of enzyme activities.However, it could also be the vector of boundary species that marks the inputs and outputs of the pathway or it could be any parameter such as a rate constant.
Obviously in practice rate constants can not be varied in a sinusoidal manner but theoretically, we could investigate what would happen if they could be varied.
The system equation ( 3) is in general nonlinear and must therefore be linearized.When linearizing, a suitable operating point must also be chosen and the most common is the steady-state.The system is at steady-state when the left-hand side of the system equation is set to zero.Therefore setting the system equation to zero gives: and linearizing using the Taylor series and truncation at the second term, yields: Note, N is considered a constant in this operation.The subscript ss is used to indicate that these values are the steady-state values.
The term ∂v(xss,p)   ∂x is a matrix of unscaled elasticities and has dimensions: r × n.We will replace this term with the less cumbersome symbol: ∂v ∂x Likewise ∂v (xss,p)   ∂p is a matrix of unscaled elasticities related to the inputs to the system.If there are p inputs then the dimensions of this matrix will be r × p.We will replace this term with the less cumbersome symbol: ∂v ∂p Using the new symbols, we can express equation ( 4) as: Those familiar with control theory will realize that equation ( 5) is a linear time-invariant system or LTI.
To examine the frequency response we must obtain the transfer function H(s).To do this we move equation ( 5) into the Laplace domain (Sauro, 2020), and assuming initial conditions are zero for the left-hand side, we solve for X(s).
Note that dδx/dt is transformed into sX(s).After rearrangement, and solving for X(s), we obtain: where I is the n × n identity matrix and s is the complex number, a + jb that the Laplace transform introduces.The equation involves a matrix inversion and we will assume that such an inversion can be done.Cases where the inverse cannot be computed include the presence of conserved moieties (which we will return to) or a model that is not realistic.
The transfer function, H(s), for a specific input and output, is defined as the ratio of the transform of the output to the transform of the input (Ogata and Yang, 2002): In terms of a matrix formulation, this would be written as: where H x (s) is the matrix form of the transfer function and includes all combinations of inputs and outputs.Comparing the above equation to ( 6) we can surmise that the transfer function must be equal to: To help clarify the notation, a subscript x has been added to the transfer function symbol H(s) to indicate that this transfer function is with respect to the species concentrations, x.
To those more familiar with control theory, this is simply the standard form obtained using the state-space approach: This form of the equation is found in many textbooks and the Wikipedia page2 offers a good summary of this equation.Matrix A is the Jacobian matrix, and comparing the standard form to equation ( 7) indicates that: A = N ∂v ∂x Before proceeding, we need to determine whether the system is stable or not.If the system is unstable then a frequency response makes no sense and no further action can be made.Stability however is easily determined by examining the eigenvalues of the Jacobian matrix.If the real part of any eigenvalue is positive then the system is unstable.
There is, however, a small complication that needs to be considered when a biochemical pathway has conserved moieties (Hofmeyr, Kacser and van der Merwe, 1986).Under these conditions, some of the eigenvalues will be zero due to row dependencies in the Jacobian.These eigenvalues should be ignored.Alternatively, some simulators, such as libroadrunner, can compute a reduced Jacobian matrix.In this case, the offending eigenvalues are automatically removed.
If a system has multiple steady-states, such as a bistable system, then there will be two sets of frequency responses, one for each steady-state.Both can be independently analyzed by picking the appropriate steady-state as the starting point.
The B term in equation ( 8) is often called the control or input matrix and is given by As indicated before, these results are fairly standard in control theory though updated to use network information such as the stoichiometry matrix and the unscaled elasticity matrix.Ingalls went further, however, and took the novel approach of using the output equation as used in control theory to describe the fluxes in transform space.In control theory, the output equation is often written as: Ingalls set the output, y, to the steady-state rates or fluxes of the pathway.We can express the reaction rate equation in vector form using: where each rate is a function of species concentrations and inputs.Linearization of this equation gives: This now looks like the output equation ( 9).Taking Laplace transforms on both sides gives: Note that the transfer function for X(s) appears in this equation.This means we can substitute equation ( 6) into this equation to give: P(s) can be moved to the right to give: The transfer function is therefore: A subscript, v, has been added to indicate that this is the transfer function for the fluxes in the pathway.
To examine the frequency response we replace the complex number s with jω3 to move us into the Fourier domain, which gives us Both equations yield complex numbers and the frequency response is encoded in these values.The phase and amplitude of the frequency response can be obtained from the corresponding phase and amplitude of the complex number.Since the frequency is given by ω, this equation can be used to determine the frequency response as a function of frequency.
Classically, the frequency response is plotted in the form of Bode plots (Ogata and Yang, 2002) which consists of two separate plots.One plot is used to indicate changes in the amplitude and a second is used to indicate changes in the phase.The x-axis in both plots marks the input frequency on a log scale which can be either expressed in radians/sec or in Hz.The amplitude plot is plotted using decibels (dB) -essentially a log scale -and the phase is plotted in degrees on a linear scale.

Generalization
Some biochemical networks, particularly signaling networks but also metabolic systems, will contain conserved moieties.This can result in numerical issues when computing the inverse because the Jacobian becomes singular.Equations ( 10) and ( 11) must be adjusted to include the link matrix (Reder, 1988), L, and the reduced stoichiometry matrix (Reder, 1988), N r .Details of this can be found in the Ingalls paper (Ingalls, 2004) although equation ( 9) in that paper has an error where the L matrix has been inadvertently left out.Note that Ingalls also only gives the frequency response equation for the independent species he calls s i .The full set can be obtained by premultiplying by L (Hofmeyr, 2001).The complete set of generalized equations is shown in ( 12).
These are the equations we will use to compute the frequency response.
Ingalls also highlighted the important fact that at zero frequency, that is H(0), equations ( 12) reduce to the unscaled control coefficients found in metabolic control analysis (Kacser et al., 1995;Heinrich and Schuster, 2012;Hofmeyr, 2001;Sauro, 2018).This result showed that MCA and classical control theory were one and the same thing though MCA derived additional theorems due to the stoichiometric nature of biochemical networks.

Example
To illustrate the use of these equations let us consider a simple two-step pathway: and compute the species frequency response using equation ( 10).We can use equation ( 10) in this example because the model does not contain any conserved moieties.
To keep things simple, we will assume that v 1 = k 1 X o and v 2 = k 2 x.This is a one-dimensional system with a single variable x but we will assume two inputs, k 1 and k 2 to make it more interesting.The differential equation for this system is: We can define the sizes for the various components using n = 1, r = 2, and p = 2. Matrix A is therefore a simple 1 × 1 matrix.The B matrix however is a 2 × 2 matrix.But as we'll see only the main diagonal has entries because k 1 has no direct effect on the second reaction and k 2 has no direct effect on the first reaction.
The A matrix is given by: Noting that ∂v 1 /∂x = 0 and ∂v 2 /∂x = k 2 we obtain: The eigenvalue for a 1 by 1 matrix is the entry itself.The eigenvalue is therefore −k 2 which is negative therefore the system is stable.
The B matrix can be computed in a similar manner, except that ∂v/∂p is a 2 × 2 matrix because there are two reactions and two possible inputs: To compute the frequency response for the system we need to insert these elements into equations ( 10) and (11).As an example, let's consider deriving H x (jω): Inserting the various terms gives: The solution has two transfer functions, one related to how changes in k 1 affect x and a second on how changes in k 2 affect x.
The frequency response can be computed by examining the magnitude and phase components of the complex numbers X o /(jω + k 2 ) and −x/(jω + k 2 ).These complex numbers need to be rationalized first into the form a + jb, from which the amplitude can be computed using the Pythagoras theorem: , Amplitude k2 = x k 2 2 + ω 2 and the phase changes using basic trig are: The phase depends on what quadrant we are in on the complex plane.The phase shift, with respect to k 1 , with vary from 0 to −90 degrees, indicating that the phase shift lags the signal.For the phase shift with respect to k 2 , there is a difference.This time the phase shift starts at +180 degrees and decreases to +90 degrees as the frequency increases.Note that the phase shift is positive for input k 2 which might suggest x is anticipating changes in k 2 .This is an illusion, because when k 2 increases, consumption of x increases so that x decreases meaning the change in x is opposite to the change in k 2 but with a delay in the response.
This example is illustrated in a later section where we demonstrate the software.
As we've seen, for small systems, it is possible to derive the analytical expressions for the frequency response but one can imagine that for large systems this can become quite tedious.The purpose of the Python package is to compute the amplitude and phase changes numerically for any sized network, with any number of inputs and outputs.

Python Implementation
Python is a good language to evaluate equations such as ( 12) because matrix multiplication and other manipulations are supported by the numpy library and Python natively supports complex numbers.As an illustration, part of the computation is shown in Listing 1.
The code begins by ensuring the system is at steady-state.If it can't find a steady-state the method throws an exception.libroadrunner has a specific method for computing steady-state as do other simulators used by the systems biology community.The variable r is a libroadrunner instance but could be another simulator that supports the features we need.
There is one minor complication when constructing the matrices.Because everything has to be computed using complex numbers we have to ensure that matrices such as the link matrix or stoichiometry matrix are converted to complex types.Moreover, libroadrunner returns named matrices so they must also be converted to pure numpy types using the numpy method array().
The evaluation of the matrix equation is fairly straightforward when using the numpy library.The code shows the individual steps as the matrix expression is constructed piece by piece.This was done to make reading the code easier and was also useful during the debugging process.
Another minor complication was computing the ∂v/∂p matrix.libroadrunner does not have a single method that can return the full matrix.The matrix was therefore constructed using a loop by calling the libroadrunner method getUnscaledParameterElasticity for each reaction.The reader may also notice that even though the unscaled parameter elasticity matrix should be a r × p matrix, the code appears to only construct a vector of size r.The reason for this is that the code as written will only compute the transfer function for a specific input.Hence we need only consider a column of ∂v/∂p that corresponds to that input.This does not preclude the possibility of looking at all possible inputs, it just means that code needs to be called multiple times with each input.Once the transfer function has been computed we use Python functions abs and math.atan2 to extract the amplitude and phase respectively.The atan2 function takes into account the quadrant to compute the correct value for the degrees.
The code is packaged into a Python function to make it easy to use.The function will compute the frequency response for a given output with respect to a given input.Two methods are provided, one for computing the species concentration frequency response and another for computing the rate or flux frequency response.
As an example, the species frequency response method has the signature: def getSpeciesFrequencyResponse(startFrequency, numberOfDecades, numberOfPoints, parameterId, variableId, useDB=True, useHz = True) The input is specified by the string argument parameterId and the output by the string argument variableId.To get the full set of transfer functions a user should call this method repeatedly with changes to these two arguments to get the full set of frequency response combinations.This does incur some inefficiency since virtually the same computation is done each time but it is often the case that only specific transfer functions are of interest.A future version might add to the code the ability to request all transfer functions to make the code more efficient.
One final aspect of the code that might be of interest is a small complication when generating the Bode phase plots.When scanning over a range of frequencies the phase shift can often go beyond -180 degrees.Normally this will result in a wraparound back to +180 degrees which when plotted appears as a discontinuity in the plot.To avoid this and to ensure a continuous appearance in the phase, the numpy.unwrapfunction is used to unwrap a signal at the discontinuous point which by default happens to be at multiples 180 degrees.Thus when the phase passes -180 degrees, the function will ensure that the phase continues beyond -180 rather than jumping back to +180.

Illustrations
With the theory and code described it is worth looking at some applications.First, we'll describe how to use the code.

Obtaining and Using the Code
The code can be obtained from https://github.com/sys-bio/frequencyResponse.There is no python pip installer because the code is just a single file called freqResponse.py.However, it is important to pip install tellurium because this is needed to supply the libroadrunner simulator and antimony support.
Once you download the freqResponse.pyfile, copy it to a location where you plan to use it.After that you need to load your model into libroadrunner.The code below illustrates this: import tellurium as te from freqResponse import * r = te.loada('''$Xo -> x; k1*Xo x -> ; k2*x k1 = 0.1; k2 = 0.34 Xo = 4 ''') Listing 2: Loading a model into libroadrunner using Tellurium where I've used the antimony syntax for describing a reaction network, but it could equally be an SBML based model.The code also imports freqResponse to make available all the methods in that file.In all the presented examples I am not going to include any explicit sink species.This is why the last reaction has an empty right-hand side.A missing sink is not important in these examples but can be included depending on the application.
The call to loada returns a roadrunner instance, r.This is required by the frequency response code.
To compute a frequency response, first, pass the libroadrunner instance to the function FrequencyResponse as follows; This will return a frequency response object, which I have called fr.Using fr one can request the frequency response with respect to species concentrations or fluxes.For example, to compute a species frequency response we call the following method: results = fr.getSpeciesFrequencyResponse(0.01, 3, 1000, 'X0', 'x3') The arguments include the starting frequency (0.01), the number of decades on the frequency x-axis (3) from the starting frequency, the number of points to generate (1000), the symbol name of the input from the model ('Xo'), and the symbol name of the output species ('x3') we wish to look at.The method for the flux frequency response has the same argument list except the output symbol must be a reaction name.In theory, I could have merged the two functions into one and the method could have recognized whether it was dealing with a species or a reaction flux.There are other optional arguments such as using Hz or rad/sec for the frequency units.
There is also a convenience method called plot which will plot the Bode plots based on the last call.The user also has access to the raw data through the results returned from the call.This will contain three columns, frequency, amplitude, and phase.

Illustrative examples
The simplest example we can try is a linear chain of reactions.It doesn't matter what the kinetic laws are because these will be automatically linearized in the algorithm.This means we can use Michaelis-Menten equations if we so wish.For illustration, however, let's just use simple mass-action kinetics.To begin we define the model using antimony (or a preexisting SBML file can also be used).A simple two-step pathway is shown below expressed as a string in antimony format:  To compute the frequency response of x with respect to input X o we use the method calls shown below.The package also provides a convenience method, plot, for plotting the results with an optional argument to generate pdf output for publication purposes.The plots illustrate two characteristics of this system, technically referred to as a first-order system because it only has one variable, x.The first is that the upper amplitude plot shows a system that behaves as a low-pass filter.That is, at low frequencies the amplitude is not changed significantly but as the frequency increases, the amplitude decreases exponentially because the system is unable to keep up with the rapid changes at the input.Since noise is usually considered high frequency, such as system could be used to smooth out a noisy signal.
The second observation is that the phase shift plot shows that x lags the input signal.At low frequencies, the phase shift is close to zero because the system can easily adjust to the changing input.This means that the input sine wave and the resulting sine wave we see at x, match each other.However, as the frequency increases the system is less able to keep up, and changes in x start to lag the input sine wave.The maximum phase change however has a limit for a first-order system.This limit is a phase shift of exactly -90 degrees.This limit to the phase shift might be counter-intuitive because one might have expected the phase to continuously shift more and more negative as the frequency increases.However, equation ( 14) confirms that the phase-shift tends to -90 10 5 10 4 10 3 10 2 10 1 10 0 10 1 10 2 0 20000 40000 Amp 10 5 10 4 10 3 10 2 10 1 10 0 10 1 10 2 Frequency 100 125 150 175 Phase (Degrees) Figure 3: The frequency response of a two-step pathway when the input is on the second step.Of note is that the phase-shift starts at +180 and decreases to +90 degrees as the frequency increases.This is because the input k 2 affects the consumption of x.
degrees at high frequencies: Why is this?Expression (14) shows that the smaller k 2 the more likely the phase shift will be -90 degrees.If k 2 is very small then we can assume that the rate of degradation is small.This means that the change in concentration of x is dominated by the input sine wave.The maximum rate of increase in x is when the input sine wave is at its maximum peak.As the input sine wave decreases the rate of increase in x slows until the input sine wave crosses its inflection point.Once the input sine wave reaches the inflection point, the level of x also peaks.Thus the input sine wave and the concentration of x will be 90 degrees out of phase with respect to the concentration of X o .This means that the delay can never be more than -90 degrees.The rate at which the phase shift approaches -90 degrees is determined by the value of the degradation constant k 2 .
In a previous section, it was also shown what the frequency response would be if the input were on the second step.The results are shown in Figure 3 which confirms what we saw with the analytical solution.
We can extend the pathway to more steps, such a model is shown below.In this case, the names of the three species are x 1 , x 2 , and x 3 .Each additional species incurs an additional maximum shift of -90 degrees.Thus x 2 could, under the right conditions, experience a maximum shift of -180 degrees.The third species, x 3 could experience an additional -90 degrees so that it could see a maximum shift of -270 degrees relative to the input X o .If you look at the right panel of Figure 4 we can see that as we move down the column, the signal becomes more and more delayed.Figure 4: Frequency responses for a four-step linear chain.The y-axis on amplitude plots has been fixed in all three cases (left column) so that the trends can be more easily seen.The amplitude drops as we move down the pathway from one species to the next.This is due to attenuation at each step.The phase plots on the right show the absolute value of the phase shift increasing by a maximum of -90 degrees across each species.This has important implications when coupled with negative feedback.

Pathway with Negative Feedback
More interesting effects can be seen when we include negative feedback.We can adjust the four-step model as shown in Listing 6.This is a linear pathway of four steps where the last species, x 3 inhibits the first reaction.This forms a simple negative feedback loop.The strength of the feedback can be changed by the Hill-like coefficient h.This system has the potential to become unstable when h is at a threshold of about 8 (Walter, 1969;Savageau, 1972).At this point, the system will destabilize and show sustained oscillations.To avoid this we will set h = 6 where the system is stable but still shows interesting frequency response behavior.Listing 6: Antimony model of a system of four steps with negative feedback from x 3 to the first step.
A frequency response of this system yields the Bode plots shown in Figure 5.The biggest change can be seen in the amplitude plot where we see a sharp spike at a frequency between 0.01 and 0.1 where.This spike is a resonance spike.
As the feedback strength gets stronger this spike becomes more pronounced until eventually the system destabilizes and switches to an oscillatory mode.This is a well-known phenomenon in control engineering and engineers try to design  systems where the resonance spike is at frequencies that the system may not experience.An alternative strategy is to reduce the strength of the feedback but this may degrade the performance of the system at lower frequencies.
Notice that the phase plot for x 3 passes through a shift of -180 degrees.This means at a shift of -180, x 3 changes in exactly the opposite direction to the input X o signal.This is the delay component that many authors some-what vaguely speak of when trying to explain the onset of oscillations in a feedback loop.
The delay itself, however, is insufficient to generate oscillations.The negative feedback loop itself will also contribute a 180 degree change in the signal.The total shift around the loop is therefore 360 degrees.This means that instead of stabilizing the system, the feedback loop coupled to the delay acts as a positive feedback loop.This is what provides the destabilizing force.If at the same time, the amplitude at the -180 crossing is one or above, the system will amplify the destabilization.This ultimately results in oscillations.A frequency response therefore gives a nice rationalization as to the origins of oscillations in a feedback loop at least for a simple feedback oscillator.Relaxation oscillators operate on a different mechanism (Sauro, 2014).

Conclusion
This article describes the theory, code, and some examples of using the software to compute the frequency response of any biochemical network.The last example in particular showed how an understanding of the frequency response in a system with negative feedback can explain the onset of oscillations.

Figure 1 :
Figure1: Injecting a sine wave into a simple two-step pathway (1) at X o .The input sine wave is shown as the blue curve, X o .The intermediate concentration of x is shown in the orange curve.Note that it takes a little time for x to stabilize to a regular sinusoidal.Once it has we can see that its amplitude is reduced and its phase is shifted.That is the peaks don't match the peaks in the blue input sine wave.The observed phase shift is due to the time it takes to fill and empty the intermediate pool x.

Figure 2 :
Figure 2: Frequency response for a two-step linear chain with respect to the intermediate x given the input X o .
import tellurium as te from freqResponse import * # load the model into libroadrunner r = r.loada(model) # pass the libroadrunner instance to the frequency response package fr = FrequencyResponse(r) # compute the species frequency response results = fr.getSpeciesFrequencyResponse(0.01, 3, 100, 'Xo', 'x') # plot the resulting Bode plots fr.plot() Listing 4: Compute the frequency response of x with respect to the input X o Running this code yields the plots shown in Figure2.

Figure 5 :
Figure5: Frequency response for a four-step linear chain with a negative feedback loop.The upper amplitude plot shows a resonance spike between a frequency of 0.01 and 0.1.
that describes the rate of change of I × r reduced stoichiometry matrix L n × n I link matrix Table 1: Matrix symbols and dimensions species in a biochemical network: